module bessels
implicit none
contains
function besseli(nu,x)
double precision nu,x,absnu
double precision besseli
double precision, parameter:: pi=3.141592653589793108624468950438d0
integer k
double precision xnu
	absnu=dabs(nu)
	k=floor(absnu)+1
	xnu=absnu-floor(absnu)

	besseli=selectk(dabs(x),xnu,k)
if ((nu < 0.0d0) .and. (xnu > 0.0d0)) then
	besseli=besseli-(2.0d0*dsin(pi*nu)/pi)*besselk(nu,x)
end if
contains
function selectk(x,xnu,k)
	double precision x,xnu
	integer k
	double precision selectk
	double precision, dimension(k)::res
	external dbsks
	call dbsis(dabs(xnu),dabs(x),k,res)
	selectk=res(k)
end function selectk
end function besseli
function besselk(nu,x)
	double precision nu,x,absnu
	double precision besselk
	integer k,NOUT
	double precision xnu
	absnu=dabs(nu)
	k=floor(absnu)+1
	xnu=absnu-floor(absnu)
	besselk=selectk(x,xnu,k)

	contains
		function selectk(x,xnu,k)
		double precision x,xnu
		integer k
		double precision selectk
		double precision, dimension(k)::res
		external dbsks
		call dbsks(dabs(xnu),dabs(x),k,res)
		
		selectk=res(k)
		end function selectk
end function besselk
function besselr(nu,x)
	double precision nu,x,ord1,ord3,dgamma
	double precision besselr
	external dgamma
	if (x .le. 1.0d-3) then
		if (nu > 0.0d0) then
			besselr=2.0d0*nu/x
		elseif (nu < -1.0d0) then
			besselr=x/(2.0d0*(-nu-1.0d0))
		else
			ord1=dabs(nu+1.0d0)
			ord3=dabs(nu)
			besselr=(dgamma(ord1)/dgamma(ord3))*(2.0d0**(ord1-ord3))*(x**(ord3-ord1))
		endif
	else
		if (dabs(nu) > 20.0d0) then
			besselr=dexp(logbesaprx(dabs(nu+1.0d0),x)-logbesaprx(dabs(nu),x))
		elseif (x > 15.0d0) then
			besselr=braket(dabs(nu+1.0d0),x)/braket(dabs(nu),x)
		else
			besselr=besselk(dabs(nu+1.0d0),x)/besselk(dabs(nu),x)
		endif
	endif
end function besselr
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function besseld(nu,x)
double precision nu,x,dgamma
double precision besseld
external dgamma
if (x < 1.0d-3) then
    if (nu > 1.0d0) then
        besseld=nu/(nu-1.0d0)
    elseif (nu < -1.0d0) then
        besseld=nu/(nu+1.0d0)
    else
        besseld=(dgamma(nu+1.0d0)*dgamma(1.0d0-nu)/(dgamma(dabs(nu))**2.0d0))*2.0d0**(2.0d0-2.0d0*dabs(nu))*exp((dabs(nu)-2.0d0)*dlog(x))
    endif
else    
    if (dabs(nu)>15.0d0) then
        besseld=dexp(logbesaprx(dabs(nu+1.0d0),x)+logbesaprx(dabs(nu-1.0d0),x)-2.0d0*logbesaprx(dabs(nu),x))
    elseif (x>15.0d0) then
        besseld=braket(dabs(nu+1.0d0),x)*braket(dabs(nu-1.0d0),x)/braket(dabs(nu),x)**2.0d0
    else
        besseld=besselk(dabs(nu+1.0d0),x)*besselk(dabs(nu-1.0d0),x)/(besselk(dabs(nu),x))**2.0d0
    endif
endif
end function besseld
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function robustlbesk(nu,x)
double precision nu,x
double precision robustlbesk

if (dabs(nu) .gt. 20.0d0) then
	robustlbesk=logbesaprx(dabs(nu),x)
elseif (x .gt. 20.0d0) then
	robustlbesk=lgbraket(dabs(nu),x)
else
	robustlbesk=dlog(besselk(dabs(nu),x))
endif
contains
function lgbraket(nu,x)
double precision nu,x,prod,j,gamr
double precision lgbraket,res,resi
external gamr
integer i
i=1
res=1.0d0
resi=1.0d0
do while (resi>1.0d-6)
	prod=.5d0+nu-i
	do j=1,2*i-1,1
	prod=prod*(.5d0+nu-i+j)
	end do
	resi=prod*gamr(i+1.0)*((2*x)**(-i))
	res=res+resi
	i=i+1
end do
lgbraket=.5d0*dlog(3.141592653589793108624468950438d0/2.0d0)-.5*dlog(x)-x+dlog(res)
end function lgbraket

end function robustlbesk
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function logbesaprx(nu,x)
	double precision nu,x
	double precision logbesaprx,root,t,argfin,pi
	pi=3.141592653589793108624468950438d0
	root=dsqrt(1.0d0+((x/nu)**2.0d0))
	t=1.0/root
	argfin=1-((3.0d0*t-5.0d0*t**3.0d0)/(24.0d0*nu))+((81.0d0*t**2.0d0-462.0d0*t**4.0+385.0d0*t**6.0d0)/(1152.0d0*nu**2.0d0))
	logbesaprx=.5d0*dlog(pi/(2.0d0*nu))-nu*root-.5d0*dlog(root)-nu*dlog(x/nu)+nu*dlog(1.0d0+root)+dlog(argfin)
end function logbesaprx
function braket(nu,x)
	double precision nu,x,dgamma,dgamr
	double precision res,resi,braket
	external dgamma
	integer i
	res=1.0d0
	resi=10.0d0
	i=1
	do while (resi > 1.0d-6)
	!print*, "h"
		resi=dgamma(.5d0+nu+i)*itgamma(.5d0+nu-i)/(dgamma(i+1.0d0)*((2.0d0*x)**i))
		!resi=tgamma(.5d0+nu+i)/(dgamma(i+1.0d0)*((2.0d0*x)**i)*tgamma(.5d0+nu-i))
		res=res+resi
		i=i+1
	end do
	braket=res
end function braket
function tgamma(x)
	double precision x,dgamma
	double precision tgamma,pi
	external dgamma
	pi=3.141592653589793108624468950438d0
	if (x>0d0) then
		tgamma=dgamma(x)
	else
		tgamma=pi/(dsin(pi*x)*dgamma(1.0d0-x))
	endif
end function tgamma
function itgamma(x)
	double precision x,dgamr,dgamma
	double precision itgamma,pi
	external dgamr,dgamma
	pi=3.141592653589793108624468950438d0
	if (x>0d0) then
		itgamma=dgamr(x)
	else
		itgamma=(dsin(pi*x)*dgamma(1.0d0-x))/pi
	endif
end function itgamma


end module bessels
